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CALCULATION  OF  GROUND  SHOCK  MOTION  PRODUCED  BY 
AIRBURST  EXPLOSIONS  USING  CAGNIARD 
ELASTIC  PROPAGATION  THEORY 


CHAPTER  1 
INTRODUCTION 


This  report  describes  a  study  which  used  elastic  wave  propagation 
theory  to  predict  and  analyze  ground  motions  produced  by  near  surface 
airburst  explosions.  The  primary  objectives  of  the  study  were  to  develop 
a  calcuiational  capability  using  the  exact  Cagniard  (1962)  elastic  formu¬ 
lation  and  to  determine  the  practical  applications  of  the  method. 

In  this  study  the  air-earth  environment  was  modeled  as  three 
homogeneous  elastic  layers — air,  soil  and  rock — separated  by  plane 
parallel  boundaries  as  illustrated  in  Figure  1.  The  air  was  treated  as 


Z 


Figure  1  Model  for  airburst  explosions  over  layered 
earth  media. 

an  elastic  fluid,  while  the  soil  and  rock  were  treated  as  elastic  solids. 
An  airburst  explosion  of  spherical  charges  was  approximated  by  a  point 
source  located  on  the  axis  of  symmetry.  Nonlinear  empirical  airblast 
arrival  time  and  overpressure  waveform  formulae  were  developed  to 


specify  the  source  characteristics  for  elastic  calculations. 


The  exact  closed  form  integral  solutions  of  Cagniard  for  the  re¬ 
flection  and  refraction  of  spherical  waves  in  elastic  solids  were 
adapted  and  extended  to  model  the  ground  shock  propagation  in  a  layered 
earth.  In  this  formulation  the  particle  motion  is  obtained  as  a  sum 
of  components  propagated  along  rays  or  paths  (such  as  shown  in 
Figure  l)  associated  with  distinct  wave  arrivals.  Calculations  using 
the  Cagniard  procedure  were  used  previously  successfully  to  predict  the 
reflection  of  underwater  explosion  shock  waves  from  the  ocean  bottom 
Rosenbaum  (1956),  Britt  (1969,  19T0),  Britt  and  Snay  (l97l)>  Snay  and 
Britt  (1973)."*"  The  theoretical  analysis  and  computer  code  development 
for  the  ground  shock  calculations  were  extensions  of  this  bottom  reflec¬ 
tion  study. 

In  the  following  sections  solutions  of  the  wave  propagation  equa¬ 
tions  for  layered  elastic  media  are  derived  using  the  Cagniard  approach. 
Much  of  the  theoretical  development  can  be  found  in  the  literature. 
Hence,  the  goal  of  this  report  is  to  present  only  the  basic  steps  of 
the  solution  procedure  and  to  bring  together  all  the  equations  needed 
for  the  ground  shock  calculations  in  a  form  tailored  to  the  problem. 

The  CAGGS  (cagniard  Ground  Shock)  computer  code  developed  for  evaluating 
these  solutions  to  obtain  particle  velocity  histories  is  discussed,  and 
comparisons  of  calculated  and  measured  waveforms  are  presented  and 
analyzed. 


Classified  reference.  Bibliographic  material  for  the  classified 
reference  will  be  furnished  to  qualified  agencies  upon  request. 
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CHAPTER  2 
THEORY 


2.1  THEORETICAL  BACKGROUND 

The  geometry  of  the  ground  shock  model  is  shown  in  Figure  1  for 
the  3-layer  case  used  in  the  CAGGS  computer  code.  Development  of  the 
theory  will  be  given  in  a  general  form  applicable  to  an  arbitrary  num¬ 
ber  of  layers.  Because  of  the  symmetry,  axisymmetric  equations  of  motion 
in  a  cylindrical  coordinate  system  are  appropriate.  The  source  is 
located  on  the  z-axis  at  the  point  (0,h)  in  the  fluid  half-space  denoted 
"air."  The  observer  (receiver  or  gage)  position  is  (r,z)  in  the  finite 
layer  of  thickness  H  denoted  "soil."  Both  the  layer  and  the  under¬ 
lying  "rock"  half-space  are  modeled  as  elastic  solids  so  that  the  labels 
"soil"  and  "rock"  are  arbitrary. 

The  solutions  of  the  elastic  wave  propagation  equations  were  ob¬ 
tained  using  the  methods  developed  by  Cagniard  (1939,  1962)  for  the 
reflection  and  refraction  of  elastic  waves  at  an  interface  between  two 
elastic  solids.  The  solutions  were  extended  to  layered  media  by  Spencer 

(i960,  1965),  Pekeris,  et  al.  (1965) ,  Abramovici  and  Alterman  (1965), 
and  Abramovici  (1970).  Additional  information  on  work  using  the  Cagniard 
approach  and  other  elastic  wave  solutions  in  layered  media  is  given  in 
the  books  by  Ewing,  Jardetzky  and  Press  (1957)  and  Brekhovskikh  (i960) 
and  the  summary  papers  of  Pao,  et  al.  (1971,  1977)  and  Britt  (1969K 
Examples  of  the  more  recent  literature  are  Abramovici  (1978),  Abramovici 
and  Gal-Ezer  (1978,  1979),  and  Pao,  et  al.  (1979). 


2.2  EQUATIONS  OF  MOTION 

Using  a  notation  similar  to  Cagniard  (1962) ,  scalar  potentials 

X  and  U  can  be  defined  for  an  elastic  solid  medium  j  such  that 
J  J 

the  radial  and  vertical  velocity  components  l  and  l  are  given  by 

r  j  zj 


(2.2.1) 
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(2.2.2) 


Z  . 
zj 


r 


These  potentials  are  to  satisfy  the  elastic  wave  propagation  equations 
in  cylindrical  coordinates 


(2.2.3) 


(2.2.4) 


where 


.  i_  +  il  ,  1  3_ 

3r2  3.2  r  Sr 


t  denotes  time,  and  c  .  and  c  .  are  the  wave  propagation  speeds  cf 

Pj  sj 

compressional  (P)  and  shear  (S)  waves,  respectively.  For  fluid  layers 

c  and  U,  are  zero.  The  subscript  j  refers  to  the  layers, 

sj  j 

The  stress  components  for  density  p.  can  be  expressed  as  follows: 

<3 


normal  stresses 


(t  )  .  =  P  . 
zz  J  J 


(Tr r>3  ’  pj 


(t,J,  =  P. 
<P<f>  J  0 


(l  -  2  c2./c2.) 
\  sj'  pj / 


2  /c2  )  +  2  c2 

sj/Cpjy/  3t2  sj  9z 


\  92X.  „  . 

!  _  2  c2  /c2  \  — i  +  2  c2  ,-P- 

c-4/Cpj)  ,t2  4  »r 


^  „  2  Vi 

+  2  c  .  — iL 
sj  r 


(2.2.5) 


(2.2.6) 


(2.2.7) 


tangential  stresses 


(t  .).  =  0 
r4>  .1 


(2.2.8) 


( 


(2.2.9) 


(t 


rz^j 


(t  )  =  0 

z4>  j 


(2.2.10) 


mean  normal  stress 


t  \  k  2  .  2 

(l  }  .  =  p  .  I  1  ~  T  c  ,/c  , 
mean  j  j  V  3  sj  pj 


at 


(2.2.11) 


In  a  fluid  the  tangential  stresses  all  vanish  and  the  normal  stresses  are 
equal  to  the  mean  normal  stress  given  by 


a2x. 


(TmQor,)f  =  -  =  P.  - O 

mean  j  j  3-^^ 


(2.2.12) 


where 


is  the  pressure. 


2.3  BOUNDARY  AND  INITIAL  CONDITIONS 


Assuming  perfect  coupling  at  the  interfaces  the  boundary  conditions 
of  the  problem  are  continuity  of  displacement  and  stress  normal  and 
tangential  to  the  interfaces.  At  an  interface  z  =  separating  media 

j  and  k  these  conditions  are 


(x 

(x 


l  .  = 

*zk 

(2.3.1) 

i  .  = 
rj 

*rk 

(2.3.2) 

zz^j 

(Tzz}k 

(2.3.3) 

rz  ^  j 

<*„>k 

(8.3.10 

In  order  that 
the  potentials 

and  in  the  limit  as 
The  source  at 


the  media  be  at  rest  before  the  source  is  initiated, 
and  U,  and  their  derivatives  must  vanish  at  t  =  0 

r  2  2*i1/2 

R  =  Ir  +  (h  -  z)  J  goes  to  infinity. 

(o,h)  is  taken  into  account  by  imposing  the 
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condition:  for  z  -*■  h  and  r  -*■  o  the  solution  tends  to  that  correspond¬ 
ing  to  a  point  pressure  (or  mean  normal  stress)  source  located  in 
medium  m 


P  = 
m 


Pm  31  (t  -  R/c  ) 


pm 


at 


--ip(t-R/c  )*  (2.3.5) 

R  o  pm 


where  X  (t  -  R/c  )  =  0  and  P  (t  -  R/c  )  =  0  for  t  <  R/c 

o  pm  o  pm  pm 

Solutions  for  other  types  of  sources  are  available  in  the  literature, 
but  the  pressure  source  is  most  appropriate  for  modeling  an  airburst 
explosion.  In  future  work  a  combination  of  various  source  types  should 
be  considered. 


2.4  PROBLEM  SOLUTION  IN  LAPLACE  TRANSFORM  DOMAIN 


In  this  chapter  we  begin  solving  the  boundary  value  problem  using 
the  Cagniard  procedure.  Briefly,  in  this  method  one  obtains  a  solution 
by  Laplace  transform  techniques.  Then  through  a  series  of  changes  of 
integration  variables  and  paths  the  solution  is  inverted  by  inspection 
back  to  the  time  domain.  The  basic  steps  and  results  are  presented 
here.  Cagniard  (1962)  and  the  other  references  can  be  consulted  for  the 
rigorous  mathematical  details. 

The  first  step  is  to  Laplace  transform  with  respect  to  time,  t  , 
the  propagation  equations  (2.2.3)  and  (2.2.1+) ,  equations  for  the 
potentials,  velocities,  stresses,  and  the  boundary  conditions.  Let  the 
transform  variable  be  s  and  denote  a  transformed  function  by  placing 
a  bar  over  the  function  symbol.  For  the  given  initial  conditions  the 
propagation  equations  become 


v2x . 


(2.1+.1) 


(2.1+.2) 


*  Note  that  the  dimensions  of  X  and  P  sire  RX  and  RF  , 

O  O  in  In 

respectively. 
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The  source  condition  (2.3.5)  becomes 
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X  (t  -  R/c  )  are  then  obtained  from 
o  pm 


-)  K 


(t  -  X)  Xjtep  (X)  dX 


(2.4.10) 


Is 

=  J  X'o  (t  -  X)  Ujtep  (X)  dX 


(2.4.11) 


Here  —  Xq  is  the  source  function.  But  in  the  equations  for  particle 

velocity  and  stress  it  will  be  advantageous  to  use  the  original  pressure 

source  function  P  (t  -  R/c  )  which  is  proportional  to  the  second 
R  o  pm 

time  derivative  of  X  in  (2.3.5). 

o 

The  significance  of  the  above  representations  is  that  solutions 
can  be  found  for  step  sources  from  which  solutions  for  an  arbitrary  time 
dependence  can  be  obtained  by  the  convolution  integrals.  The  solutions 
Xjtep  and  U^tep  will  be  treated  in  the  following  sections.  For 
simplicity  of  notation  the  superscript  "step"  will  be  omitted,  except 
when  necessary  for  clarity,  and  it  will  be  understood  that  the  convolu¬ 
tions  must  be  performed  for  a  specified  source  function. 

Transformed  solutions  in  any  of  the  layers  can  be  written  in  the 


00 

'  /* pj(u 


)  exp(sot  .z)  J  (sur)du 
pj  o 


QO 

/vu) 


(2.4.12) 


+  /  b  .(u)  exp(-sa  .z)  J  (sur)du 

/  pj  pj  ° 


00 

=  Jasj(u 


)  exp(sct  .z)  J'(sur)du 
sj  o 


(2.4.13) 


+  /  b  ,(u)  exp(-sa  ,z)  J'(sur)du 

/  sj  sj  o 


where 


/  2  -2\1/2  /  2  -2\1/2 
°P.3  \U  CPJ/  ’  asj  "  V  sjj 


2  -2  \ 

+  c  I  ,  and  J  ( sur )  i s 
SJ  /  o 

the  Bessel  function  of  the  first  kind  of  order  zero.  The  functions  a 


and  b  are  to  be  determined  from  the  boundary  conditions.  The  terms 
containing  a^  and  a^  represent  P  and  S  waves,  respectively, 
traveling  in  the  negative  z  direction.  Similarly,  the  b  and  bgj 
terms  correspond  to  waves  traveling  in  the  positive  direction. 

The  transform  of  the  step  source  can  be  put  in  a  similar  form 
using  the  Sommerfeld  integral  derived  in  Ewing,  et  al.  (1957) 


iR  exp(-sR/cpm)  =  f JQ(sur)  exp(-sa  |h  -  z|)  ^ 
J  Pm 


(2.U.1U) 


The  direct  solution  for  the  coefficients  a  and  b  is  very  com¬ 
plicated  even  for  the  3-layer  case  and  results  in  very  involved  formulae. 
(See  Abramovici  and  Alterman  (1965)  for  example.)  An  alternative  ap¬ 
proach  is  the  generalized  ray  path  concept  introduced  by  Spencer  (i960, 

1965).  It  was  shown  that  the  potentials  X  and  U  can  be  built  up 

J  J 

from  terms  which  represent  distinct  arrivals  or  rays  in  the  form 


=  X  XJq  =  J  jV)  exp^-s(Sokdk)J  Jo(sur)du 


(2.U.15) 


VI  Ujq  =  ^ /Cq(u)  expj^-s(Eak  V]  Jq( sur )du 


(2.U.16) 


where  q  is  an  identifier  of  the  ray,  Cq(u)  contains  factors  represent¬ 
ing  the  source  and  the  generalized  reflection  and  transmission  coeffi¬ 
cients  for  each  interface  encounter.  The  subscript  k  is  used  symbol¬ 
ically  to  identify  segments  of  the  ray  q  ,  and  d,  is  the  vertical 
projection  of  the  ku  segment.  If  a  segment  is  in  medium  m  ,  then 
denotes  and  asm  for  P  and  S  waves,  respectively.  A  ray 

represents  a  term  in  X.  or  U  if  the  last  segment  which  joins  to  the 

J  J 

receiver  is  a  P  wave  or  S  wave,  respectively. 


Spencer  (i960)  derived  the  method  of  generalized  rays  by  first  con¬ 
sidering  the  interactions  at  a  single  interface  for  incident  P  and  S 
waves.  These  interactions  can  be  identified  with  two  letters:  the 
first  denoting  the  type  (P  or  S)  of  the  incident  wave  and  the  second 
the  type  of  the  reflected  or  transmitted  wave.  There  are  then  four  com¬ 
binations  PP  ,  PS  ,  SP  and  SS  for  both  reflections  and  transmissions 
for  a  total  of  eight  possible  interface  interactions  for  which  coeffi¬ 
cients  must  be  determined.  Once  derived  the  formulae  can  be  used  at  any 
interface  by  inserting  the  applicable  wave  speeds  and  densities. 

Starting  at  the  source,  rays  can  be  drawn  representing  all  the 
possible  types  of  paths  which  join  the  source  to  the  receiver.  The 
appropriate  coefficients  (derived  in  the  next  section)  are  then  included 
in  C^(u)  for  each  interface  encounter  of  the  ray.  By  considering  the 
ray  segment  beginning  at  the  source  which  has  no  interface  contact, 
the  "source  factor"  (derived  in  Section  2.6)  to  be  included  in  C  (u) 

q 

can  be  determined. 

The  arrival  time  (derived  in  Section  2.1l)  along  each  ray  depends 
on  the  wave  speeds  and  the  number  and  lengths  of  segments.  At  a  given 
time  only  those  rays  which  "have  arrived"  must  be  considered. 


2.5  DERIVATION  OF  GENERALIZED  REFLECTION 
AND  TRANSMISSION  COEFFICIENTS 

In  this  section  equations  are  derived  for  the  generalized  reflec¬ 
tion  and  transmission  coefficients.  Consider  the  case  of  waves  incident 
only  from  medium  j  onto  medium  k  .  Let  the  interface  separating  media 
j  and  k  be  z  =  0  .  Substitute  the  solutions  in  the  form  of  (2.4.12) 
and  (2.4.13)  into  transformed  versions  of  the  boundary  conditions  (2.3.1) 
-(2.3.4)  using  (2.2.1),  (2.2.2),  (2. 2. 5)-(2. 2. 10)  to  express  displace¬ 
ment  and  stress  components  in  terras  of  the  potentials.  Assuming  that 
medium  j  is  above  k  ,  the  terms  containing  a^^  and  a^ ^  represent 
the  incident  waves  and  hence  b  ^  =  bg^  =  0  .  The  resulting  equations 
can  be  put  in  a  form  similar  to  that  used  by  Ewing,  et  al .  (1957)  for 
plane  wave  reflection  using  the  notation 


'  ck  c=k/"3c'll 


(2.5.1) 


for 


i  =  j  or  k 


(2.5.2) 


2  -2 

Ai  =  ^  +  Csi 

h  -  («*k  -  *?)  ^-5.3) 


^  *  JTJ  '  V  'sj 

(2.5.U) 

L  =  -fc  (a,  -  2Qu2)  c2. 

3  “pj  V  j  /  SJ 

(2.5.5) 

H  ’  2“V  (1  -  «>  %i 

(2.5.6) 

Mi  "  2“V  <3  -  «>  '1 

(2.5.7) 

-  sj  (AJ  -  ^ )  'l) 

(2.5.8) 

=  <«**- V  csj 

(2.5.9) 

Mu  *  b 

(2.5.10) 

After  some  manipulation  the  boundary  conditions  yield  equations  as 
follows : 


L,  a  ,  +  M,  a  ,  =  a  .  +  b  . 
1  pk  1  sk  pj  pj 


Vpk  *  "sS*  ■  asJ  -  \s 

(2.5.12) 

L3apk  +  M3ask  "  apJ  "  bpj 

(2.5.13) 

Llapk  *  V.k  ■  asj  +  \j 

(2.5.1**) 

lU 


Setting  =  0  ,  agj  is  the  only  incident  term  and  equations 
(2.5.H)-(2.5.lM  can  be  solved  for  the  generalized  reflection  coeffi¬ 
cients  K  and  transmission  coefficients  T  for  an  incident  S  wave. 
Define 


D  =  (1^  +  L3)(M2  +  M^)  -  (L2  +  LU)(M1  +  M  )  (2.5.15) 


to  obtain 


KJ3H  =  ^=2(L3M1-Ll“3,/D 

SJ 


(2.5.16) 


k”  -  -  [(L2  -  LU)(M1  +  M3)  -  +  L3)(M2  -  Mu)]/D  (2.5.1?) 


T1k  =  ^  =  "2(M1  + 


(2.5.18) 


(2.5.19) 


where  the  superscripts  denote  the  types  of  waves  before  and  after  the 
interface  interaction,  and  the  subscripts  denote  waves  incident  from 
medium  J  onto  medium  k  . 

Incident  P  waves  can  be  considered  by  setting  agj  =  0  .  The 
resulting  coefficients  are 


=  iEi=  [(Lj  -  L3)(M2  +  M^)  -  (L2  +  Ll4)(M1  -  M3)]/D  (2.5.20) 


(jk  “  ~  2{HM2  "  L2M1i)/D 

r  J 


(2.5.21) 


(2.5.22) 


. -Et.  .  2(„2  .  Ml|)/D 


Tjk  =  r7'-2<L2*H>/D 

P  J 


(2.5.23) 


Since  both  media  j  and  k  are  solids,  coefficients  and 

T^j  for  waves  incident  from  media  k  can  be  obtained  by  simply  inter¬ 
changing  the  subscripts  in  the  above  equations.  In  the  case  that  one 
medium  is  a  fluid  and  one  is  a  solid  this  symmetry  does  not  hold,  but 
the  coefficients  can  be  obtained  by  taking  the  limit  as  cg  of  the 
fluid  goes  to  zero. 

Denote  the  fluid  medium  by  j  and  the  solid  k  .  After  taking 

the  limit  c  -*■  0  ,  the  coefficients  above  can  be  simplified  by 
s  J 

defining 


f  ■  "A 


Djk  *  aPj(4  -  ‘m%kask)  *  f°; 


pkcsk 


(2.5.2U) 


(2.5.25) 


The  resulting  coefficients  for  waves  incident  from  the  fluid  are 


K3k  ’  [°pj(4  -  '•“"V’sk)  - 


f0pkcsk  ^DJk 


(2.5.26) 


^k  ~  2fapjAkcsk^DJk 


(2.5.27) 


=  Uufa  ,a  ,  c  ?/D, . 
jk  pj  pk  sk'  Jk 


(2.5.28) 


For  waves  incident  from  the  solid  medium  k  the  coefficients  are 
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Kkj  =  '‘‘"“pjWV 

=  -[“pj(4  +  bvv)  *  fapkcsk]/BJk 

KJ’! =  -[aPj(\2  *  ka\k%i)  -  faPkcsk]/Djk 

*w  “  -'"“pjvvv 

'IkJ  =  2apkAkcsk//Djk 


(2.5.29) 

(2.5.30) 

(2.5.31) 

(2.5.32) 

(2.5.33) 

(2. 5.3*0 


The  coefficients  for  shear  components  in  the  fluid  are  all  zero. 

For  the  case  of  two  fluid  media  the  coefficients  simplify  to 


K??  =  (a 


jk 


Pj 


-  fV)/uPj 


+  fa  .  ) 

pk 


(2.5.35) 


2fa 


PJ 


fa  ,  ) 
Pk' 


(2.5.36) 


Again  there  is  symmetry  for  waves  incident  from  medium  k  ,  and  now  all 
coefficients  involving  shear  waves  are  zero. 


2.6  DERIVATION  OF  THE  SOURCE  FACTOR 


The  source  factor  in  C  (u) 

q 


can  be  determined  by  equating  the 


i 
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integral  representation  (2.4.14)  for  a  P  wave  source  with  the  q  =  0 
(no  interface  encounters)  term  of  (2.4.15) 


00 

/  Co(u 


)  exp(-sa  d  )  J  (sur)du 
pm  o  o 


(2.6.1) 


=  /J o(sur)  exp("sapmih  '  Z')  SHf 


Equating  expressions  inside  the  integrals  gives 

C  exp(-sa  d  )  =  — —  exp(-sa  |h  -  z|  ) 
o  v  pm  o  V  pm1  '/ 

Let 


(2.6.2) 


pm 


d  = 
o 


|h  -  zj 


(2.6.3) 


and 


I 


C  =  u/a  (2.6.4) 

o  pm 

Thus  C  (u)  for  a  P  wave  source  can  he  expressed  in  the  form 
<1 

C  (u)  =  —  H  (KT)  (2.6.5) 

pm 

where  II  (KT)  denotes  a  product  of  reflection  and  transmission  coeffi- 

q 

cients  for  each  interface  encounter  of  ray  q  . 

2.7  A  RAY  NUMBERING  SCHEME  FOR  3-LAYER  MEDIA 

The  3-layer  geometry  of  Figure  1  was  used  for  the  ground  shock 
calculations  in  this  report.  A  simple  ray  numbering  scheme  for  this 
case  is  as  follows.  First  group  together  rays  having  the  same  number  of 
reflections  within  the  layer.  Then  within  each  group  arrange  rays 


)■ 

I 
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t 


according  to  an  increasing  number  of  S  wave  segments.  S  segments 
extending  the  full  thickness  of  the  layer  are  counted  higher  than  seg¬ 
ments  which  cross  only  part  of  the  layer.  This  scheme  arranges  rays 
within  a  group  according  to  increasing  arrival  times  for  the  common 

situation  in  which  |z|  <  H/2  and  c  <  c  /2  . 

1  s  p 

Each  ray  can  be  identified  by  the  letters  P  or  S  of  its  seg¬ 
ments.  The  first  segment  in  the  fluids  is  always  P  and  its  letter 
designation  will  be  dropped.  The  last  segment  crosses  only  part  of  the 
layer  and  hence  has  a  valve  of  d^  ,  vertical  projection,  which  is  less 
than  H  .  The  other  segments  in  the  layer  have  =  H  .  For  calcula- 
tional  efficiency  rays  having  the  same  Ea^d^  (and  hence  the  same  arri¬ 
val  times)  are  numbered  together  and  computed  with  the  same  integration 
since  all  terms  in  the  integrals  are  the  same  except  for  the  product  of 
reflection  and  transmission  coefficients  HKT  in  C  (u)  .  Table  1 

q 

lists  rays  1  through  12  with  the  factors  Ea^d^  and  the  products  ITKT  . 
Note  that  rays  9  and  10  are  degenerate,  having  two  components. 

For  higher  ordered  reflections  the  labels  PnSmS  or  PnSmP  can 
be  used,  where  n  and  m  denote  the  number  of  P  and  S  segments, 
respectively,  crossing  the  layer.  The  last  letter  denotes  the  type  of 
the  last  segment.  Using  this  notation  the  general  formula  for  Ea^d^ 
for  ray  q  is 


lW  =  hapl  +  nHap2  +  ^82 


(2.7.1) 


The  choice  of  the  upper  or  lower  terms  in  braces  depends  on  the  last  seg¬ 
ment  of  the  ray.  H  +  z  is  used  if  the  segment  is  upward  (n  +  m  odd), 
and  -z  is  used  if  the  segment  is  downward  (n  +  m  even),  is  used 

if  the  last  segment  is  a  P  wave  (q  odd),  and  a ^  used  for  a  S 

wave  (q  even). 

2.8  EXAMPLE  OF  THE  CONSTRUCTION  OF  TRANSFORMED  SOLUTIONS 


As  an  example  of  how  the  ray  solutions  X,  and  U 


are 


Table  1.  Characteristic  parameters  of  rays  q.  =  1  to  12 


1 


Path  diagram 

Path  type 

nATk 

1 

P 

a  _  h  - 

-  a 

tpp 

\ 

pi 

p2 

12 

* 

2 

S 

a  .  h  - 

-  a  z 

T?* 

Pi 

s2 

12 

3 

h 

5 

6 


PP 

PS 

SP 

ss 


aplh  +  “p2(2H  +  z) 
aplh  +  ap2H  +  “s2(H  +  z) 
aplh  +  “s2H  +  ap2(H  +  Z) 
aplh  +  “s2(2H  +  z) 


TPP^PP 

T12K23 

mPPv-PS 

T12K23 

TPsKsp 

T12K23 

mPSj.SS 

T12K23 


t.1 


7 

8 
9 

10 


PPP 

PPS 

PSP 

SPP 

PSS 

SPS 


vh  +  V(2H " z) 

aPlh  +  “p2(2H)  '  “s2Z 
aplh  +  ap2CH  "  Z)  +  as2H 


a1h+aH+a0(H-z) 
pi  p2  s2 


11 

SSP 

aplh  "  ap2Z  +  as2(2H) 

12 

sss 

aplH  +  as2(2H  "  z) 

mPPt/PPi/'PP 
X12K23  21 

mPPKPP^PS 
12  23  21 

mPPKPS„Sp 

T12K23K21 

Tps  sp  PP 
T12K23K21 

mPPlfPS^SS 

T12K23K21 

ps  sp  ps 
T12K23K21 

TPSICSSICSP 

T12K23K21 

„PS  SS  SD 
T12K23K21 
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constructed,  consider  the  path  B  shown  in  Figure  1  which  is  denoted 

q  =  5  in  Table  1.  Since  the  last  segment  is  a  P  wave  in  medium  2, 

this  ray  represents  the  term  X„  .  The  ray  is  first  transmitted  from 

ps 

medium  1  into  medium  2  with  a  coefficient  T£2  .  The  ray  is  then  re¬ 
flected  at  the  interface  z  =  -H  with  a  coefficient  •  Hence,  from 


23 


(2.6.5)  we  obtain 


p  -  _JL_  tPs  wsP 
°5  ~  a  12  K23 
Pi 


(2.8.1) 


The  vertical  projections  dL^  of  the  ray  are  h  ,  H  and  H  +  z  , 
so  that 


£  aA  =  aplh  +  as2H  +  V(H  +  Z 


(2.8.2) 


the  solution  for  the  ray  is  then 


-5"/[v^K”]'wH0plh 


+  a^H  +  ap2(H  +  Z)J  f  JQ( sur )du 


*1 


(2.8.3) 


Solutions  for  the  other  rays  are  built  up  in  the  same  manner . 

2.9  INVERSION  OF  THE  TRANSFORMED  SOLUTIONS 

General  terms  of  (2.U.15)  and  (2.U.16)  are 


03 

"jq  =/ e*P  [-s(Zokdk)]  Je(sur)du 


(2.9.1) 


00 

p  q=J  Cq(u)  exp  [-s(lakdk)]  j;(sur)du 


(2.9.2) 


These  expressions  can 


be  inverted  using  Cagniardrs  method  to  obtain  the 
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solution  X.  and  U.  in  the  time  domain.  The  first  step  is  to 

jq  jq 

substitute  into  (2.9-1)  and  (2.9.2)  the  integral  representations 


J  (sur)  =  Re 
o 


,/2 

'I 


exp(-isur  cos  io)dw 


(2.9.3) 


J'  (sur) 
o 


=  Im 


«/2 

H 


cos  a)exp(-isur  cos  w)dw 


(2.9.U) 


where  Re  and  Im  denote  real  and  imaginary  parts,  respectively,  and 
i  =  /-l  .  Rearrange  the  equations  to  obtain 


jq 


7h  "f  r 

=  J  \  ^  Re  J  C^(u)  exp  -s^£  +  i'ur  COS 


doAdu  (2.9.5) 


jq 


°r(  */-2  p 

=  J  Im  J  C^(u)  cos  a)  exp  -s  d“A  +  iur  cos  uj 


du>  du  (2.9.6) 


The  next  step  is  the  key  Cagniard  substitution:  a  change  of  integration 
variable  from  u  to  t  at  constant  u>  defined  by 

(2.9.7) 


t  = 


4.  •  1 

to  give 


^k&k  +  iur  cos 


v/2 


jq 


=  j~<~Re  j  Cq^u^(|t’) 


(2.9-8) 


Several  rigorous  mathematical  details  have  been  omitted  here.  See 
the  references  (Cagniard  (1962)  or  Britt  (1969)  for  example)  for  addi¬ 
tional  information. 
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cos  a)  dcu>exp(~st  )dt 


(2.9.9) 


^=/{  -1"/ vu)(£i 

From  the  definition  of  the  Laplace  transform  the  expressions  in  braces 
are  the  inverse  transforms 


X.  =  -  Re 

jq  it 


(2.9.10) 


tt/2 

U.  =  -  Im  /  C  (u)(!t}  cos  a)  dm  (2.9.11) 

jq  it  I  q  W 

Jo 


This  completes  the  essential  steps  of  the  potential  solutions.  The 
following  expresses  the  equations  in  forms  more  suitable  for  computa¬ 
tions  . 


2.10  CHANGE  OF  INTEGRATION  VARIABLE  FROM  w  TO  u 

The  fact  that  u  in  the  integrands  of  (2.9.10)  and  (2.9.11)  can¬ 
not  in  general  be  explicitly  determined  as  a  function  of  m  makes  the 
numerical  integration  of  these  equations  difficult.  A  change  of  inte¬ 
gration  variable  from  w  to  u  at  constant  t  produces  expressions 
more  amenable  to  computations 


(2.10.1) 


EVk)du 


(2.10.2) 
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where  we  have  defined 


u2r2  +  (t  -  %  ak^ 


(2.10.3) 


and  ,  corresponding  to  u>  =  0  ,  is  the  solution  in  the  fourth 
quadrant  of  the  complex  plane  having  the  smallest  modulus  of 


t  =  Iak(u^)dk  +  iu^r 


(2.10.U) 


and  u^  corresponding  to  u>  =  ir/2  is  the  solution  of 


t  =  rak(u2)dk 


(2.10.5) 


The  imaginary  parts  of  integrals  (2.10.1)  and  (2.10.2)  are 
symmetric  about  the  real  axis  and  the  real  parts  are  anti-symmetric. 
The  integrals  from  u^  to  u^  can  be  replaced  by  integrals  from  u^ 
to  u^  where  u^  is  the  complex  conjugate  of  u^  .  The  results  are 


X,  =  ~ 

jq  lir 


1  f  Cq(u)du 


(2.10.6) 


u-, 

iirr  I 


)(t  -  Eakdk) 


(2.10.T) 


2.11  ARRIVAL  TIMES 

An  examination  of  the  integrands  above  reveals  that  the  integrals 
are  zero  before  a  time  tn  which  will  be  determined  later.  In  addition 

iq 

there  is  a  characteristic  time  t^  >  t^  which  is  the  maximum  value  of  t 

for  which  u^  from  (2.10.1*)  is  a  pure  imaginary  number  u^  =  -iy^  for 

v  >  0  .  This  time  t„  corresponds  to  the  arrival  time  of  a  ray 

J/2  2q 

satisfying  Snell's  Law  of  Refraction.  The  condition  necessary  for  the 


maximum  is  obtained  by  applying  3t/3y2  =  0  for  Up  =  -iy  in  (2.10.!+): 

-  r/r  .  0  (2.11.1) 


This  gives 


*2»  *  E\(-iyanM>  dk  *  y2max  r 


(2.11.2) 


We  can  identify 


y2max  =  Ck1  sin  6k 


(2.11.3) 


where  0^  is  the  angle  the  k  segment  of  the  ray  makes  with  the 
normal  to  the  interface.  This  equation  is  really  a  statement  of  Snell's 
law  since  k  ranges  over  all  segments  of  the  ray  and  y  is  the 


same  for  all  k  . 


We  next  examine  the  case  in  which  t  <  as  defined  by  (2.11.2). 

An  inspection  of  the  integrands  (2.10.6)  and  (2.10.7)  shows  that  the  integrals 
are  zero  until  at  least  one  of  the  square  roots  becomes  imaginary 
or  complex.  It  can  be  shown  that  u^  and  u^  are  zeros  of  y  which 


means  the  potentials  are  non-zero  for  t  >  t„ 

-  2q 


For  t,  <  t0 
lq  2q 


exist,  one  of  the  ex's  in  the  integrands  must  have  a  zero  at  a  value 

of  u  =  +iy,  such  that  y,  <  y_  .  But  it  has  been  shown  that  y- 
—  1  1  2max  J  2n 

is  always  less  than  any  of  the  c  which  occur  in  la  <1  .  Hence, 

k  K  k 

/  t  o  c- 


y^  must  be  the  zero  of  some 


/  -1  2V 

a  =  C  _  y  I 
pn  \  pn  ^l/ 


occurring  in  C^(u) 


but  not  in  la,  d,  .  (Since  c  <  c  ,  only  a  need  be  considered.) 
k  k  pn  sn  J  pn 


Denote  the  minimum  such  c  ^  as  c  ^ . 

pn  pmm 


Thus  t,  <  t„  occurs  if 
lq  2q 


’  .  <  y„ 

pmin  2max 


u,  =  -ic  .  ,  equation  (2.10.!*)  gives 

1  pmm  a 

t,  =  la,  (ic  1.  )  d,  +  c"1.  r 

lq  k\  pmm/  k  pmm 


(2.11.U) 


In  the  following  discussion  we  generalize  this  definition  of  t,  to  be 

-1  4 

the  minimum  arrival  time  of  ray  q  ,  hence  if  no  cpmin  ^2max  exlctr'’ 

then  we  set  t,  =  t0  . 

lq  2q 


2.12  THE  CHANGE  OF  INTEGRATION  VARIABLE  FROM  u  TO  v 


An  additional  change  of  integration  variables  in  the  potential 
solutions  (2.10.6)  and  (2.10.7)  is  advantageous  because  y(u)  is  zero 


at  the  integration  limits  u^  and  -Q.^ 


It  is  easily  shown  that  the 


-1/2 


and 


integrands  have  singularities  which  behave  like  (u  -  u..  ) 

-1/2 

(u^  -  u)  ,  except  at  t  =  when  a  logarithmic  singularity  occurs. 

The  first  step  in  eliminating  the  singularity  for  t  ^  tg  is  to 
replace  the  integration  from  u^  to  u^  by  one  in  the  first  quadrant 
of  the  complex  plane.  Using  a  path  of  constant  real  part,  Re(u^)  , 
the  calculation  can  be  performed  using  a  real  integration  variable  if 
we  set  u  =  Re(u^)  +  iy  to  give 


(2.12.1) 


(2.12.2) 


where  y„  =  Im(un )  and  y  =  c  \  for  t  <  t_  and  y  -  0  for 
2  X  1  pmin  cq  ± 

t  >  *2q  ' 

Next  make  a  change  of  variable  (similar  to  that  used  by  Longman 
(1961))  from  y  to 


or 


w 


(2.12.3) 


Re(u)  =  y  =  y2 


2 


w 


(2.12.IO 
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to  give 


(2.12.5) 


(2.12.6) 


where 


w1  =  (y2  -  y1) 


1/2 


Since  the  quotient  w/y  has  a  finite  limit 


as  w->*0(y->-y),  the  singularity  has  been  eliminated. 


2.13  SOLUTIONS  FOR  VELOCITY  COMPONENTS 


In  explosion  effects  research  particle  velocity  (or  velocity  ob¬ 
tained  by  integrating  acceleration)  and  stress  are  usually  the  quanti¬ 
ties  measured.  Equations  for  velocity  and  stress  components  are  best 
derived  by  performing  the  various  space  and  time  derivatives  on  the 
transformed  potential  solutions  (2. 9-1)  and  (2.9.2)  rather  than  deter¬ 
mining  the  derivatives  after  inverting  these  solutions. 

Since  velocity  v  is  related  to  the  displacement  £  by 
v  =  d£/dt  ,  the  transform  is  v  =  st  for  zero  intial  displacement. 

Equations  for  the  transformed  velocity  components  can  then  be  obtained 

from  the  transformed  displacements  by  simply  multiplying  by  s  .  Simi- 

2  2  2— 
larly,  the  3  X,  /3t  terms  in  the  stresses  have  transforms  s 

jq  Jq 

Denote  the  P  wave  contribution  to  the  radial  velocity  for  ray  q 

s 


Similarly,  for  the  vertical 


as  v13  and  the  shear  wave  term  v‘ 

rq  n  s  rq 

velocity  use  tt  and  v  .  Let  S  =  +1  for  rays  in  which  the  last 
zq  zq  q 

segment  is  downward  and  S  =  -1  if  the  last  segment  is  upward.  Velocity 
components  are  positive  upward  and  outward.  Substitute  the  transformed 
potentials  (2.9.1)  and  (2.9.2)  into  transformed  versions  of  equations 
(2.2.1)  and  (2.2.2),  and  multiply  by  s  to  obtain 


— p 

vr  - 
rq 


3X, 

= 

3r 


OO 

.*/  u  C^(u)  expf-sla^d^]  <V(sur)du 


(2.13.1) 
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V 

rq 


3U, 

_ Ja  = 

3z 


oo 

7%o 


.S  C  (u)  exp[-sZa,  d,  ] 
J  q  q  k  k 


J' ( sur )du 
o 


(2.13.2) 


/ 


zq 


ax, 

3z 


a  .S  C  (u)  exp[-sEa,  d,  ]  J  (sur)du 
pj  q  q  k  k  o' 


(2.13.3) 


o 


Note  that  these  equations  are  solutions  for  a  step  wave  source  from 

which  solutions  for  an  arbitrary  potential  time  variation  XQ(t  -  R/c  ) 

can  be  obtained  by  using  (2.k.l0)  and  (2.U.11). 

2 

The  s  factors  above  can  be  eliminated  if  a  pressure  source 

P  =  7  P  (t  -  R/c  )  is  used  instead  of  the  potential  source.  By 
m  R  o  pm  _  _ 

applying  equation  (2.k.3)  which  relates  Xq  and  Pq  ,  it  is  easily  seen 
for  a  velocity  component  v  that  we  can  define 

v  =  -  P  V  (2.13.5) 

P  o 
m 

where  V  is  obtained  from  equations  (2.13.l)-(2. 13.^ )  by  dropping  the 

p 

s  factor  up  front.  Then  v  is  obtained  from  the  convolution 


v 


t 

/P’(t  -  X)  V(X)  dX 
° 

o 


(2.13.6) 


Note  that  subscripts  and  superscripts  have  been  dropped  for  simplicity. 

The  V  equations  involve  J’(sur)  as  in  U.  and  can  be  in- 
r  o  _  jq 

verted  in  the  same  manner.  Similarly,  the  equations  contain  Jo(sur) 

and  can  be  inverted  as  X,  was.  The  results  analogous  to  (2.12.5)  and 

jq 

(2.12.6)  are 
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where  S  =  +1  or 
<1 

upward,  respectively. 


-1  when  the  last  segment  of  ray  q  is  downward  or 


2.1U  SOLUTIONS  FOR  STRESS  COMPONENTS 

For  a  pressure  source  we  can  obtain  step  wave  response  functions 
Z  and  Y  for  the  stress  components  which  are  analogous  to  the  response 
V  for  the  velocity  components.  Then  corresponding  to  (2.13.6)  the  stress 
components  can  be  obtained  from 


t 

-J-|  p- 

J  ° 


(t  -  X)  Z  (X)  dX  - 

Pm 


t 

/ 


r  (t 

o 


X)  Y ( X )  dX  (2. 1)4.1] 


where  the  proper  superscripts  and  subscripts  are  implied  on  both  sides 
of  the  equation.  The  following  stress  response  functions  were  derived 
using  the  same  procedure  used  for  the  V  functions: 


4 


29 


I 


7sq  _  k  2 

L  -  -  p  .  c 

rz  Trr  j  sj 


wi  r 

f  4 


A  C  (u) 

■* -  (t  -  Za  cL  )  wdw 

uy  k  k 


(2.1U.H) 


zpq 

mean 


= *  "jI1  - 1  c=j,cpj)  /  Re[-4^] 


(2.  .12) 


All  Y  response  functions  are  zero  except  Y1 4  ,  Y  4  ,  Y,4  ,  and  Y^4 

rr  rr  <p<J>  <j>4> 

which  are  listed  above-  In  addition,  7- SC1  ,  Z  ??  ,  and  7 sq  are  zero. 

v  »  AA  *  '  m  nn  ki 


2.15  LIMITING  VALUES  OF  THE  POTENTIAL  SOLUTIONS  AT  t  =  t„  AND  r  =  0 

2q 

The  case  for  which  c  ^ .  <  y~  =  c.-1  sin  0,  corresponds  to  the 

pmin  2max  k  k  - 

supercritical  reflection  of  plane  wave  theory.  This  case  occurs  when 

sin  0,  >  sin  0  . ^  =  c,  /e  .  .  It  can  be  shown  for  the  supercritical 

k  cnt  k  pmin 

reflection  that  the  responce  to  a  step  wave  source  has  a  logarithmic 

singularity  at  t  =  t^  •  However,  if  the  source  function  and  its  time 

derivative  vanish  more  rapidly  than  a  logarithm  diverges  as  t  -►  F/c  , 

the  response  obtained  after  convolution  is  finite. 

For  subcritical  cases  (sin  0.  <  sin  0  ..  )  the  responses  at  t  =  t„ 

k  cnt  c(\ 

have  finite  values.  The  stepwave  response  functions  for  the  potentials 
(2.12.5)  and  (2.12.6),  for  the  velocities  (2.13.7)-(2.13.10) ,  and  for 
the  stresses  (2.lL.2)-(2.lL.12)  can  all  be  written  in  the  form 


F(t)  = 


/ 


(2.15.1: 


Then  using  the  procedure  of  Cagniard  (3962)  pages  bfi  and  89,  the  value 


of  F  at  t  =  t. 


for  r  >  0  is 


*1  k  v 


F(t_  )  = 


Tm[f(u,t  )] 

_ gq 


rj  -2.3 
y„  r  Xd,  c.  / a 
2max  k  k  k 


tl/2 


(2.15-2) 


evaluated  at  u  =  iy_  where  y__  and  t„  are  riven  by  (2.11.1) 
J2max  2max  2q 

and  (2.11.2),  respectively. 

In  the  case  r  =  0  ,  the  solutions  simplify  considerably.  The 
expressions  for  the  potentials  can  be  easily  derived  by  setting  r  =  0 
in  equations  (2.9-1)  and  (2.9-2)  and  using  ,T  (0)  =  1  and  .'P(O)  =  0  . 
Equation  (2.9.7)  reduces  to 

t  =  Zci^d^  (2.15.3) 


and  u  is  real  rather  than  a  complex  variable.  Then  making  the  change 
of  integration  variables  from  u  to  t  as  in  section  2.9,  the  solutions 
can  be  inverted  to  obtain 


X.  (r  =  0) 

Jd 


u£V\ 


(2.15.1*) 


U.  (r  =  0)  =  0  (2.15.5) 

jq 


where  u  is  obtained  from  (2.15.3).  The  arrival  time  arguments  of 
Section  ? .11  do  not  apply  for  r  =  0  .  For  this  case  the  arrival  time 
is 


t,  =  t0  =  (2.15.6) 

lq  2q  k  k 

which  results  from  u  =  0  .  Tn  (2.15.1*)  the  limit  an  u  -*■  0  of  C  (u)/u 
must  he  taken  at  the  arrival  time. 

The  corresponding  velocity  and  stress  stepwave  responses  for  a 
?  wave  source  are  of  the  form 


F(r  =  0) 


(2.15.7) 


j  f(u) 

uWk/o,k 


where  the  iiotation  of  (2.15.1)  has  hec*n  used.  Note  that  all  response 

terms  containing  t  -  in  the  numerator  vanish  at  r  =  0  since  this 

factor  results  from  .T'(sur)  which  is  zero  at  r  =  0  .  All  other  velocity 

o 

and  stress  component  terns  involving  2  waves  are  also  zero  at  r  =  0  since 
they  result  from  U.  and  its  derivatives. 
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CHAPTER  3 


THE  CAGG3  CODE 


3.1  GENERAL  DESCRIPTION 

A  computer  program  CAGGo  (. Cagniard  ground  Chock)  was  written  to 
calculate  the  particle  velocity  components  v^  and  v^  for  the  three 
layer  ( fluid-solid-solid)  case  where  the  point  source  is  in  the  fluid 
and  the  receiver  is  in  the  1'inite  layer.  Displacement  calculations  ob¬ 
tained  by  numerically  integrating  the  velocities  are  also  available.3' 

The  code  can  calculate  up  to  183  arrival  times  which  include  all 
rays  having  12  or  fewer  reflections  and  the  P33P  ray  having  13  reflec¬ 
tions  .  However,  reflection  and  transmission  coefficient  products  II KT 
are  coded  for  only  71  of  the  typically  earliest  arrivals.  The  available 


rays  are 
p^p  -  p^sp 


p  -  i  °! 


p°p 


P3S3S 


pfp  -  r 


AOy,  p  9„(, 

r  -  r  o  o  , 


pup  -  p10p>r 


5s2p 

12, 


PbP  -  fVs 


-  pi;lss 


13 

and  P  P  . 


These  ray.-,  allow  cal  •ulating  velocity-time  curve.:  for  approximately  3  to 
U  2  wave  transit  time:',  of  the  layer  after  the  first  arrival. 

The  CAGGP  code  war.  written  in  iuRTRAN  for  the  Honeywell  0OO-6OOO 
series  computer.; .  The  code  is  run  in  the  time  sharing  mode  from  a 
Tektronix  **0i '■*  high  s:  ecu  graphics  terminal.  Hard  copy  of  the  informa¬ 
tion  on  the  screen  is  available  at  the  push  of  a  buttom.  The  output  is 
a  tab  Le  V'  ferity  an  i  d : .:  placement  versus,  time  and  velocity  and  dis- 
history  plots  such  as  the  calculated  curves  in  Figures  3-10. 
Tie  curve.-  r.iii  wn  hav*'  beer;  reduced  to  about  1/2  the  original  size. 


.a.  input; ■  .-.No  ■jtput: 


'able  2  .-.hows  an  examj  Le  calculation  including  input  and  tabular 
Oatnut  evry  tent:  1  time  'i’he  decay  exponent  allows  calculations 

for  a  1/r  ‘  decay  in  layer  .  .  AIL  computations  in  this  report  used  a 
pure  linear  decay  with  an  exponent  1.0.  The  next  input  selects  which 
plots  art  les  :  re.: .  .hr  gu  *«»  frecueucy  input  allows  calculation  of 


A  variation  ,>f  t.ne  ''AG  >lc  which  oalculnteu  stress  components  in 
medium  in  al  .  available. 


1 


3l> 


d 


motion  record:’,  which  r. xmulato  the  rcr.ponuo  of  velocity  rare:;  to  the 
theoretical  velocity.  The  variable  J.W  specifies  the  type  of  explosive 
device.  IW  -  1  for  nuclear  and  I W  =  ?  for  conventional  explosives 
(H.E.  ) .  W  is  the  charge  weight  in  ki.Lotons  nuclear  or  pounds  T.N.T. 
equivalent  weight  li.E.  The  next  two  lines  are  the  wave  speeds  {ft/msec) 
and  densitites  (gr:i/em^)  in  media  2  and  3.  the  parameters  in  the  air 
(medium  i)  are  not  input  but  are  determined  by  the  program  from  empiri¬ 
cal  formulae.  NM1N  and  WRAY  specify  the  minimum  and  maximum  values 
of  q  to  be  used.  H  ,  h  ,  r  ,  X  are,  respectively,  the  layer 
thickness,  explosive  source  height  of  burst,  horizontal  range,  and  -z 
the  depth  of  gage  or  receiver.  These  variables  are  specified  in  feet. 
The  next  two  inputs  are  the  site  altitude  and  air  leg;  erature  used  in 


,  t  1  on .  i  lie 


ay  where 


calculating  the  airblast. 

The  first,  line  of  output  gives  the  arrival  ?  and  peak 

overpressure  (nsi)  at  z  =  0  directly  a  novo  t.l-  a  •  i  t  i  on.  The 
following  table  gi  ves  characteristic  parameters  f  ■  •.  ray  whore 

IR  =  q 

T1  =  time  of  arrival  (msec) 

1  1 

T2  =  time  of  arrival  (msec)  t , 

''-1 

TCONST  =  time  constant  G.,  (msec)  described  later.  For  the 
present  version  of  CAGGS  the  e  |_t}  are  the  same  for 
all  q  . 

RW  =  horizontal  distance  r  at  which  the  ray  enters 
medium  2  ^ 

PAIR  =  airblast  peak  overpressure  (psi)  at  range  KW 

DEFAULT  DT  =  the  time  step  (msec)  which  is  adequate  for  most  runs. 

With  experience  larger  steps  can  often  be  used. 

TSTOP  =  1 .  ime  (msec)  of  the  end  of  the  calculation.  The  code 
presently  allows  up  to  1000  steps  and  quits  automati¬ 
cally  if  the  maximum  is  reached. 

The  next  output  is  a.  table  containing  the  stop  mun.ber  IP  ,  time  t 
(msec),  velocities  v  and  v  (ft/sec),  and  displacements  P  ^  and 
?r  (ft).  The  maximum  and  minimum  values  of  the  velocities,  displace¬ 
ments,  and  the  velocities,  measured  by  a  simulated  velocity  gage  are 


age  are 


listed  following  the  table. 


All  of  the  ’nput— output  doner . bed  above  is  printed  on  the  terminal 


screen.  Hard  copy  is  obtained  by  pressing  the  "copy"  key.  If  plots 
were  specified,  the  user  presses  the  "return"  key  after  the  copying  is 
completed.  The  program  clears  the  screen  and  draws  a  plot.  The  code 
stops  after  each  plot  to  allow  copying  and  continues  when  the  "return" 
key  is  depressed.  After  a  case  is  completed  the  programs  allows  chang¬ 
ing  all  or  only  part  of  the  input. 


?. 3  NUMERICAL  INTEGRATION  PROCEDURE 


To  obtain  particle  velocity  components  and  ,  two  numeri¬ 

cal  integrations  must  be  performed  for  each  ray.  First,  the  stepwave 
responses  must  be  evaluated  from  equations  (2. 13. T )-( 2. 13. 10 ) .  Second, 

the  convolution  integrals  in  the  form  of  (2.13.6)  must  be  computed. 

To  start  the  calculation,  y.  is  computed  for  each  ray  from 

2max 

(2.11.1)  using  the  Method  of  False  Position  (Conte  (1965))  combined 

with  an  increment  halving  iteration  scheme  to  insure  convergence  in 

difficult  cases.  Arrival  times  t,  and  t^  are  then  determined  from 

lq  ^q 

(2.11.2)  and  (2.11 .b)  as  appropriate.  The  source  amplitude  and  other 

time  independent  terms  of  the  step  responses  are  calculated  in  the  same 

section  of  the  code.  The  stepwave  responses  are  then  evaluated  using  a 

four  point  Gaussian  quadrature  with  eight  subintervals  for  a  total  of 

32  points.  The  and  integrals  are  calculated  together  since 

most  of  the  factors  in  the  integrals  are  the  same.  The  integration 

limit  Q  is  obtained  from  a  Newton's  Method  iteration  (Conte  (1965)) 

for  u,  in  (2.10.U).  When  t,  <  t  <  t„  ,  the  substitution  u,  =  -iy„ 
1  lq  2q  12 

is  used  to  allow  calculation  with  real  arithmetic.  When  t  >  t0  , 
complex  arithmetic  must  be  used  in  the  iteration.  In  either  case,  an 
initial  guess  based  on  the  value  of  from  the  previous  time  step 

generally  produces  rapid  convergence. 

The  convolution  integration  is  simplified  by  expressing  the  source 
pressure  time  dependence  P  (t)  in  the  form 


-U-b/c  )/ek 

a.  e  1 


(3.3.1) 


37 


x+Ax  x  Ax 

Because  of  the  property  e  =  e  e 


the  convolution  can  be  written 


t 

/I"  it  -  A)  V  (  A  )  d  A 

° 

o 


t-At 


-(t-At-A)/0 

V(A)dA 


(3.3.2) 


+ 


) 

t-At 


-(t-A  )/0R 


V  ( A  )d.\ 


if  P  =  H ( t  -  H/c  )  p(t  -  R/c  )  where  P(0)  ^  0  ,  we  must  add  an 

additional  term  to  (2.13.6)  correspotidinc  to  ■  ll(t.  -  h/Cp^)  ■ 

6(t  -  R/e  ,  ).  1'lie  result  is 
pi 


o 


t^(t  -  A)  V( A )dA  =  p (0)  V(t)  + 


o 


p'(t  -  A)  V ( A )dA 


(3.3.3) 


The  last  integral  can  then  be  expressed  as  in  (3-3.2). 

The  first  integral  in  (3-3.2)  in  the  convolution  from  the  previous 
time  step  multiplied  by  an  exponential.  Thus  in  a  time  step  we  must 
only  calculate?  the  second  integral.  For  times  not  near  a  singularity 
time  t,,,  for  a  supetvr  if.  ical  reflection,  we  integrate  an  follows. 
Denote  t,j  =  t  -  At  nr.: 

V1  =  V(t  -  At)  (3.3-^) 

V,  =  V(t)  (3.3.5) 

Use  linear  int.<?rpo  L.ati.x.  i  et.ween  the  calculated  points. 


(3.3.6) 


V  =  V1  +  (V2  -  Vi)(t  ~  V/At 


then 


I 


2 


/ 

t-At 


-(t-X)/0 

e  V(X)dX 


(3.3.7) 


Near  a  singularity  time  t0  we  make  the  change  of  variable 

x  =  ^2q  "  X!1/2  (3.3.8) 

to  eliminate  the  logarithmic  singularity  from  the  integral  producing 


+2 


2  - 


(t) 

c 

-(t-x)/e 

I 

e  k  V( X  ) 

L  J 

xdx 


x( t-At ) 


(3.3.9) 


where  +  and  -  apply  to  t  >  tn^  and  t  <  ,  respectively.  This 

integral  is  evaluated  using  a  generalized  Simpson's  rule  obtained  by 
integrating  a  parabola  fit  to  V(t)  ,  V(t  -  At )  ,  and  V(t  )  where 
t  is  the  time  of  the  step  before  t  -  At  .  A  fixed  time  step  is  not 
required.  The  convolution  integrations  (3-3.7)  or  (3.3*9)  take  only  a 
small  percentage  of  the  computing  time  required  to  evaluate  the  V(t)  , 
hence  several  terms  in  the  sum  of  (3.3.1)  can  be  used  to  fit  the  source 
time  dependence  Pq ( t )  without  adding  significantly  to  the  total  com¬ 
puting  time. 

3.1*  MODELING  THE  AIRBLAST  WAVEFORM 

The  point  source  amplitude  and  waveform  used  in  the  CAGGS  code  wan 
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chosen  to  simulate  the  airblast  overpressure  produced  by  airburst  explo¬ 
sion  of  spherical  charges.  Empirical  curve  fits  based  on  explosive 

field  test  data  were  used.  Peak  overpressure  1’  and  arrival  time 

rn«.xx 

t  at  the  air-soil  interface  (z  =  0)  for  height  of  burst  h  and  hori- 
a 

zontal  range  r  were  obtained  from  formulae  extracted  from  the  ANSWER 
(Analysis  System  for  Weapons  Effects  Research)  code  (Britt  (1980) ) .  The 
nuclear  section  uses  the  Brode  (197D,  1978)  height  of  burst  curves.  The 
conventional  explosive  (H.E.)  values  are  based  on  curves  developed  for 
the  ANSWER  code. 

For  the  ranges  of  interest  for  most  linear  ground  shock  calcula¬ 
tions  and  for  the  common  surface  tangent  burst  configuration  (charge 
resting  on  the  ground),  the  H.E.  peak  overpressure  and  arrival  time  curves 
can  be  approximated  by 

P  -  3770  r  psi  for  2.5  <  r  <  1?  ft /lb  (3.^.1) 

t  ~  O.OUk  r1'9-5  msec/lb1/3  for  2  <  r  <  15  ft/lb1/3  (3.4.2) 

where  W  is  the  explosive  charge  weight  in  pounds  T.N.T.,  r  =  r/W1^3 

and  t  =  t  /W3^3  .  The  value  of  P  given  here  is  that  measured  at 
a  a  max 

ground  level  and  hence  contains  both  the  incident  (or  source  contribu¬ 
tion)  and  the  reflection.  Because  of  the  great  impedance  mismatch  be¬ 
tween  the  air  and  soil,  the  linear  theory  predicts  a  reflection  virtually 
identical  to  the  incident  pulse  if  the  soil  layer  is  very  thick.  In 

order  to  match  the  total  pressure  (normal  stress  1  )  at  the  interface, 

zz 

one  half  of  the  empirical  blast  overpressure  is  vised  for  the  incident 

pressure  in  the  linear  calculations. 

The  airblast  waveforms  P  (t)  such  as  shown  in  Figure  2  were  fit 

o 

to  a  sum  of  exponential  terms  of  equation  (3.3.1)  using  data  from  the 
Department  of  Defense  (I'.O.D.)  explosive  tests  series  given  in  table  3- 
The  events  used  in  the  waveform  model  were  mostly  surface  tangent  bursts 
supplemented  b,y  other  near  surface  detonations.  The  parameters  a  and 

K 

0^  were  chosen  t.o  fit  the  experimental  parameters  noted  on  Figure  2  as 
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Figure  2  Typical  airblast  waveform. 
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Table  3.  Explosive  test  series  used  in  the  airblast  model. 


Test  Series 

Explosive  Charge 

CENSE  1 

1000  lbs  NMa 

CENSE  2 

300  lbs  TNT 

CENSE  3 

200  lbs  NM 

Dial  Pack 

^00  tons  TNT 

Mine  Ore 

100  tons  TNT 

Mine  Under 

100  tons  TNT 

Middle  Gust  3 

100  tons  TNT 

Middle  Gust  4 

100  tons  TNT 

Mineral  Rock 

100  tons  TNT 

Misers  Bluff  Phase  I 

1000  lbs  TNT 

Misers  Bluff  Phase  il 

120  tons  ANFO^ 

Mixed  Company 

500  tons,  20  tons  TNT 

Pre-Dice  Throw  11 

100  tons  TNT,  120  tons  ANFO 

Pre-Mine  Throw  IV 

1004  lbs  TNT,  '(.2  tons,  102.  k  tons  NM 

1969  Height-of-Burct  Series 

(BRL) 

1000  lbs  TNT 

NOTE :  A  table  f  factors  for  converting  inch-pound  units  of  measurement 
metric  fO  i)  units  is  presented  -n  page  3- 

Nitromethane  (NM) 

Ammonium  nitrate  and  fuel  oil  (ANFO) 


It? 


well  as  the  positive  phase  impulse  I  (the  integral  of  the  curve  from 
t  to  t  ) .  A  large  body  of  data  was  obtained  from  reports  on  P  , 
1^  ,  t  and  t+  .  A  representative,  but  much  smaller,  quantity  of 
information  on  the  negative  phase  was  used.  More  data  was  available 
from  the  agencies  conducting  the  tests  or  from  D.O.D.  archives,  but  time 
and  funding  limitations  prohibited  the  extensive  use  of  this  material. 
The  scatter  of  the  data  used  in  determining  the  negative  phase  param¬ 
eters  is  much  greater  than  for  the  positive  phase.  It  is  expected  that 
the  negative  phase  model  used  here  will  not  be  changed  greatly  when  a 
larger  data  base  is  used. 

The  resulting  waveform  for  T.N.T.  explosions  is 


P  (t)  =  P 
o  max 


(3-1*.  3) 


where  P  ^  is  the  peak  overpressure  in  psi ,  W  is  the  charge  weight 

in  pounds  of  T.N.T.  ,  t  =  t/W1/3  msee/lb17"3  ,  0  =  0  /W^3  msec/lb1/3  , 

K  K 


a0  =  minimum 


.>00  P 


0.1*8 

max 


{  1*010 


(3.14.1*) 


al*  = 


[21  P0'1  for  5  <  P  <  29-1  psi 
1  max  ~  max  — 


7930  p-1-66  for  P  >29.1  psi 
max  max 


(3.1*.  5) 


/  6.26  P  por  5  <  p 


or  >  <  l’  <  200  pr*i 
max  — 


6,  =  <218  P  1,3t>  for  200  <  P  <  700  psi 
1  I  max  ^  max  — 


0.081*  P_0'lb  for  P  >  7 00  psi 
max  max 


U.*4.6) 


P  feiiWBSSW&Sv-S** 


1;3 


II  P  U,‘L°  for  9  <  P  <  200  psi 
max  ~  max  — 


6,  =  \  ?kkb  P-0'13  for  200  <  P  <  TOO  psi 
k  I  max  max 


(3.h.  ■() 


O.U9  for  P  >  TOO  psi 
L  max 


0g  =  minimum 


0.5  9, 


(3.^.8) 


0.3  =  0.99  B0 


(3.U.9) 


0,,  =  l.Oi  0, 

9  h 


(3.^ .10) 


A  finite  rise  time  can  be  introduced  by  adding  an  additional  term 


to  (3.U.3)  of  the  form  a^  e 


,  where 


ar  =  -a. 
u  1 


(3-U.ll) 


Let  At  be  the  scaled  rise  time  of  the  pulse,  then  if  the  last  four 
o 

terms  of  (3.1*>3)  change  only  slightly  in  this  time,  0^  can  be  approxi- 
mat ed  by 


(3. *4. 12) 
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The  amplitude  P  of  P  (t)  is  then  obtained  bv  replacing  a,  =  1’ 

max  o  1  max 


-At/ 0  -At/ Qr 
1  6 
e  -  e 


(  3 .  !• .  1 3 ) 


3.5  APPROXIMATION  OF  THE  AIR BLAST  WITH  A  POINT  SOURCE 


As  can  be  seen  from  the  empirical  formulae  described  above,  the 

airblast  cannot  be  directly  modeled  as  a  linear  point  source  in  our  range 

of  interest  (10  <  P  <  200  psi)  since  (a)  the  airblast  peak  overpres- 
~  max  — 

sure  decays  like  l/R^-dT  instead  of  1/R  ,  (b)  the  blast  propagation 
rate  decreases  with  range  instead  of  the  constant  speed  c  ,  and 
(c)  the  shape  of  the  waveform  changes  with  range  rather  than  being  a 
fixed  function  PQ(t  -  R/ Cpj_ )  •  Within  the  framework  of  the  Cagniard 
model  two  options  are  available  for  approximating  the  empirical  airblast: 
(a)  simulating  the  pressure  in  space  and  time  by  a  distribution  of  point 
sources  or  (b)  using  a  single  point  source  chosen  to  linearize  around 
a  particular  range  related  to  the  time  of  dominant  motion.  The  first 
option  would  be  more  accurate  but  would  likely  require  much  more  com¬ 
puter  time  to  evaluate.  The  simpler  and  less  costly  linearization 
approach  (b)  was  used  in  this  study.  Several  procedures  were  investi¬ 
gated,  but  linearization  around  ray  q  =  2  ,  the  directly  transmitted 
shear  wave  (path  A  of  Figure  l),  produced  the  best  overall  agreement  with 
measured  velocity  waveforms  for  materials  ranging  from  weak  soils  to 
hard  rock. 

The  directly  transmitted  P  wave  path  was  determined  by  iterating 
for  the  range  rg  at  which  the  ray  enters  the  soil.  The  empirical 
formulae  (such  as  {?,.b.2))  were  used  to  calculate  the  blast  arrival  time 


t  fo.  an  initial  value  of  r  =  r 
a  s 

in  the  air  was  computed  from 


An  average  1’  wave  nneed  c 


(■■  •  -r 


(3.5.1) 


Then  from  equations  (2.11.2)  and  (2.11.3)  we  obtained  the  relation.'. 


■  ('1  -  £  )1/2 "  *  (■-;  -  /■:  )' 
\  pi  2max J  \  s,  .  max/ 


+  y.. 


r*l  *  s 


^  2max 


(3.1).  3 


which  must  also  he  satisfied  for  the  ray.  Those  equations  are  actually 

a  disguised  statement  of  Snell's  Law  of  He fraction.  Suing  the  initial 

values  of  r  and  e  .  in  (3.5.T)  and  (i.[q.3)  a  new  estimate  of  r 
s  ul  s 

was  obtained.  It  is  easily  seen  that  the  solution  lies  between  r  and 
this  second  estimate  of  ro  .  Hence,  iteration  was  continued  using  an 
interval  halving  technique  on  this  interval  until  (  3.  5  •  1 )-(  3  ■ ‘i .  3  )  were 
satisfied  within  an  acceptable  tolerance.  The  point  source  amplitude 
and  pulse  shape  were  then  chosen  using  the  empirical  formulae  to  match 
tiie  airblast  at  the  point  (r  ,0).  The  resulting  c  and  point  source 

5  pi 

function  P  (t  -  R/c  )/R  were  then  used  for  all  the  rays, 
o  pi 


CHAPTER  4 


COMi’AK  I  t'.OK  OF  CALCULATED  AND  MEASURED  PARTICLE 
VELOCITY  WAVEFORM:' 

4.1  SELECTION  OF  TEST  OKR1EO 

Calculations  were  performed  for  the  three  CENSE  (Coupling  Effi¬ 
ciency  of  Near  Surface  Explosions)  explosive  field  test  series  ( Ingram 
1.07Y,  1080).  These  tests  were  chosen  because  they  provide  a  variety 

of  site  characteristics  which  were  relatively  well  controlled.  The 
test  beds  were  either  effectively  homogeneous  or  had  layering  suitable 
for  the  two-layer  model.  In  audition,  eacii  of  the  series  had  near  sur¬ 
face  airburst  explosions  for  which  particle  velocity  or  acceleration 
was  measured  in  the  upper  layer  for  a  variety  of  ranges  from  the 
explosions . 

4.2  CALCUliATl ONU  FOR  CENSE  1 

CFNC.K  1  consisted  of  a  series  of  1000-lb  spheres  of  nitromethane 
( ll40  lbs  TNT  equivalent)  detonated  over  a  massive  Kayenta  sandstone 
formation.  These  events  provide  data  for  chocking  the  calculations  for 
motion  in  a  strong,  homogeneous  material  which  behaves  elastically  for 
stress  levels  of  hundreds  of  psi.  The  surface  rock  was  thick  enough 
that  three  layers  wore  not  needed  for  the  computations.  Figure  3  com¬ 
pares  the  theoretical  and  measured  vertical  and  radial  velocity  compo¬ 
nents.  Vertical  velocity  is  positive  for  upward  motion  and  radial 
velocity  is  positive  for  outward  motion.  Note  that  the  experimental 
and  theoretical  curves  arc  plotted  on  different  scales  and  that  the 
calculation  represents  only  part  of  the  measured  curve.  The  material 
properties  used  for  this  calculation  were  p  =  0.0112  gm/cra  , 

e  0  =  9  ft /msec  ,  e  ,  =  4  ft  ./msec  and  p.  =2.4  gm/cm  Event.  1. 

P1-  t: 

shown  in  the  figure,  vn  :  detonated  with  its  charge  center  0  ft  above 

the  rock.  Measurement  r.  were  made  with  velocity  gages  (having  a  nominal 

600-Hz  frequency  response)  placed  2  ft  U  low  the  rock  surface.  The 

airblast  peak  pressure  above  the  gage  was  measured  at  ~>1  psi  at  the 


r 
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W  ft  range.  Figure  4  shows  a  calculation  at  the  36-ft  range  (120  psi). 
Only  the  horizontal  velocity  record  was  obtained.  There  is  good  agree¬ 
ment  between  calculations  and  experiment  at  both  ranges  except  in  the 
magnitude  of  the  second  positive  phases.  This  slight  discrepancy  is 
probably  a  result  of  the  simple  source  linearization  and  not  an  effect 
of  non-linear  material  properties  of  the  rock. 

4.3  CALCULATIONS  FOR  CENSE  2 

The  data  of  CENSE  2  provided  information  to  check  the  elastic  cal¬ 
culations  for  a  two-layered  clayey-silt  soil  site.  The  explosives  used 
were  300-lb  spherical  TNT  charges.  In  the  calculations  shown  in 
Figures  5  and  6,  the  soil  was  modeled  in  two  iayers:  a  surface  layer 
20  ft  thick  with  c^2  =1.1  ft /msec  ,  cg2  =  0.6  ft /msec  ,  and  p2  = 

1.7  gm/enf5  ,  and  a  lower  half-space  with  c  _  =  1.6  ft /msec  ,  c  = 

,  p3  s  3 

0.7  ft/msec  ,  and  =  1.75  gm/cnr  .  The  measurements  were  made  with 
velocity  gages  located  1.5  ft  below  the  surface.  Event  2  vas  detonated 
with  charge  center  7.2  ft  above  the  soil  surface.  Excellent  agreement 
was  obtained  in  Figure  5  at  the  67- ft  (13-jisi)  range.  Figure  6  at  the 
43-ft  (34-psi)  range  shows  slightly  poorer  agreement  but  still  within 
typical  scatter  of  field  measurements.  At  a  range  of  32  ft  (6 0  psi), 
Figure  7,uhe  linear  calculations  begin  to  fail  to  reproduce  the  major 
characteristics  of  the  measured  motion.  In  this  case  the  linear  theory 
does  not  predict  the  large  initial  downward  and  outward  displacements 
(area  under  the  velocity  curve)  seen  it:  'he  experiments.  These  differ¬ 
ences  are  probably  a  result  of  the  nonlinear  material  properties  of  the 
soil  becoming  important  and  a  result  f  *  no  close-in  source  conditions 
not  being  adequately  modeled  by  tin  1  eon  1  i  -el  airilas.t.  input,  used  in  the 
calculations. 


4.4  CALCULATIONS  FOR  CENSE  3 

CENSE  3  provided  measurements  for  eompuri.  '  .  with  t  n<  cry  !\ 


weak  soil  layer  over  a  hard  rock  side.  This  series,  eons  i  si  <‘d  of  s<-v<  i: 
explosions  of  200  lb  (226  lb  TNT  equivalent )  of  n  i  treret  bane .  T*,.-  * 

bed  consisted  of  compacted  backfill  of  "alluvium"  soil  t inm  d  ■  v  r  a 


Figure  4  Comparison  of  particle  velocity  calculations 
with  CENSE  1  measurements  at  36-ft  range. 


CENSE-?  EVEN’ 


ii-ure  7  Comparison  of  particle  velocity  calculations 
with  CEIiSE  2  measurements  at  32- ft  range. 


The  thickness  of 


Kayenta  sandstone  deposit,  similar  to  that  of  CENSE  1. 

tiie  soil  was  varied  from  0  to  u  ft.  Measurements  of  vertical  and  radial 

acceleration  were  made  at  middepth  in  the  soil  layers  and  in  the  rook. 

Velocity  histories  were  obtained  by  integrating  the  acceleration  records 

Events  2  and  k  were  surface  tangent  bursts,  that  is,  tiie  explosive 

charge  was  resting,  on  the  soil  surface.  The  soil  layer  thickness  in 

Figure  8  was  6  ft.  In  Figure  9  the  thickness  was  3  ft.  Tiie  material 

properties  used  in  the  calculations  were  c  „  =  0.9  ft/moec  c  _  = 

P2  s2 

0.3  ft/msec,  and  p„  =  1.6  gm/em  for  the  soil  and  c  _  =  8  ft /msec  , 

2  p3 

c  _  =  3  ft/msec,  and  p  }  =  2.h  gm/ern  for  the  sandstone. 

At  the  56-th.  (12-psi)  range  of  Figure  8  tiie  calculations  are  in 
good  agreement  with  tiie  experimental  curves  up  to  a  time  of  about 
h 5  msec  if  the  high  frequency  spikes  are  neglected.  These  spikes  re¬ 
sult  from  using  an  airblast  pulse  with  zero  rise  time  and  reflections 
at  a  perfectly  distinct  interface.  High  frequency  motion  of  this  type 
is  filtered  out  of  tiie  measurements  because  of  tiie  nonlinear  effects 
of  tiie  soil  and  finite  frequency  response  if  tiie  gages  and  recording 
system.  At  the  32- ft  (52-psi)  range  of  Figure  9  agreement  is  poorer 
but  tiie  initial  velocity  amplitudes  are  still  close  to  tiie  measured 
values.  At  a  range  of  ft  (IT)  psi)  not,  shown  the  calculated  initial 

peaks  are  nearly  a  factor  or  two  higher  than  the  experimental. 

Figure  10  is  presented  for  comparison  with  Figure  9  to  illustrate 
the  effect  of  further  increasing,  tiie  soil  layer  thickness  at  the  CENSE  3 
site.  The  initial  portions  of  the  records  are  produced  by  the  directly 
transmitted  ['  and  S  waves  and  are  not  iependent  on  tiie  soil  thick¬ 
ness.  The  later  motion  is  a  e.  .mpl i cat od  interaction  of  reflected  waves 
for  moderate  layer  thickness.  in  going  from  a  layer  thickness  of  3  ft 
as  in  Figure  8  to  tiie  6- ft  lay-  in  Figure  d  the  change  in  frequency 
of  the  motion  is  roughly  prop  ed  ional  to  tiie  layer  thickness  change, 
but  tiie  waveforms  at  12- ft  thickness  do  not  follow  this  pattern. 

From  the  few  waveforms  presetted  here,  one  cannot  draw  general 
conclusions  on  factors,  affecting  t lie  period  and  amplitude  of  the  low 
frequency  motion.  A  detailed  parameter  study  and  analysis  will  be  neces 
sary  to  determine  how  the  motion  ciiang.es  in  going  from  very  thin  to  very 
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thick  layers.  It  appears  that  simple  rules  of  thumb  based  on  S  or 
p  wave  layer  transit  times  will  be  valid  in  only  very  restricted 
ranges  of  thickness  and  elastic  parameters. 
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CHAPTER 


CONCLUSIONS  AND  RECOMMENDATIONS 

Calculations  performed  with  the  CAGGS  code  were  in  good  agreement 
with  measurements  from  the  three  CENSE  test  series  which  had  very  differ¬ 
ent  site  conditions  and  particle  velocity  waveforms.  These  sites  span 
much  of  the  spectrum  of  typical  soil-rock  environments.  CENSE  1  site 
was  a  hard,  strong,  thick  sandstone  rock;  CENSE  2  site  was  a  two 
layered  soil  medium;  CENSE  3  site  was  a  layer  of  weak  soil  backfill  over 
a  thick  sandstone  formation.  Based  on  these  calculations,  we  conclude 
that  predictions  of  velocity  waveforms  using  the  Cagniard  formulation 
of  the  elastic  theory  and  the  localized  airblast  source  model  can  be 
expected  to  be  accurate  within  the  scatter  of  explosive  tests  measure¬ 
ments  for  times  up  to  about  one  cycle  of  the  low  frequency  motion  and  for 
airblast  overpressure  levels  at  the  gage  range  up  to  approximately  Uo 
psi  for  explosions  over  weak  soils  and  over  100  psi  for  strong  rocks. 

At  pressure  levels  in  the  range  of  hO  to  100  psi  for  explosions  over 
soil,  the  elastic  theory  still  predicts  the  general  character  of  the 
motion  but  overestimates  the  peak  velocities  and  underestimates  the 
large  initial  downward  and  outward  displacements. 

The  simple  localized  airblast  source  model  linearized  around  the 
directly  transmitted  shear  wave  may  be  a  major  contribution  to  failure 
of  the  calculations  at  higher  pressures  and  at  late  time  on  waveforms. 
Future  work  should  be  directed  at  improving  the  source  formulation, 
but  any  model  other  than  a  single  point  source  would  greatly  increase 
the  computer  time  costs  of  the  calculations. 

The  linear  wave  propagation  model  can  produce  motion  waveforms 
for  homogeneous  sites  such  as  the  CENSE  1  sandstone  for  very  small  com¬ 
puter  costs.  But  because  of  the  rapidly  increasing  effort  required  to 
calculate  the  late  time  portion  of  waveforms  for  layered  media,  routine 
calculations  to  times  greater  than  about  three  shear  wave  transit  times 
of  the  layer  appear  to  be  more  expensive  than  comparable  linear  finite 
difference  methods  or  normal  modes  techniques.  Calculations  with  the 
Cagniard  theory  for  more  than  two  soil  layers  appear  to  be  quite 
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expensive  except  for  early  time  motion  or  special  cases  where  reflections 
in  one  layer  can  be  neglected  or  modeled  by  a  few  rays.  Our  incomplete 


calculations  have  indicated  that  three  solid-layer  cases  such  as  present 
at  the  Dial  Pack  site  where  a  shallow  water  table  exists  with  a  much 
deeper  rock  layer  can  be  effectively  modeled  by  the  Cagniard  method 
using  only  a  few  carefully  selected  rays  in  the  thin  surface  layer. 

The  primary  applications  of  the  CAGGS  code  for  computing  motion 
waveforms  in  layered  media  appear  to  be  early  time  motions  up  to  about 
two  shear  wave  transit  times.  These  cases  require  relatively  low  com¬ 
puting  costs  and  minimal  effort  to  change  code  inputs.  Runs  to  later 
timer;  and  for  more  than  three  layers  are  practical  but  are  considerably 
m  're  expensive.  Since  the  theory  follows  rays,  the  composite  waveforms 
•can  le  dissected  to  study  the  contributions  of  individual  arrivals. 

“'his  property  of  the  method  makes  It  ideal  for  studying,  the  basic  char¬ 
acteristics  and  effects  of  the  controlling  parameters  of  wave  propaga- 
ti  >nc  in  layered  media.  The  CAGGS  code  developed  in  this  study  appears 
to  be  an  ideal  to'.l  for  performing  detailed  parametric  studies  of  ground 
shock  at  intermediate  ranges  where  the  airblast  pressure  is  from  5  to  50 
psi  in  weak  soils  anu  to  hundreds  of  psi  for  strong  rock. 
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